# Replication Material for "Social Capital, Institutional Rules 
#     and Constitutional Amendment Rates"
# William Blake, Joseph Cozza, Dave Armstrong and Amanda Friesen
#
# before running, execute setup.r first
#
# This file produces tables Tables A4 and A5 and Figures 1, A4 and A5. 

## Cross-national Models
cndat <- import("crossnational_replication.dta")
set.seed(207)
m1n <- brm(bf(amendments ~ vp + words + govConfid + offset(log(age)), 
              shape ~ vp + govConfid), 
           data=cndat, 
           family=negbinomial(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars = save_pars(all = TRUE))


m2n <- brm(bf(amendments ~ vp + words + partyConfid + offset(log(age)), 
              shape ~ vp + partyConfid), 
           data=cndat, 
           family=negbinomial(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars = save_pars(all = TRUE))
m3n <- brm(bf(amendments ~ vp + words + courtConfid + offset(log(age)), 
              shape ~ vp + courtConfid), 
           data=cndat, 
           family=negbinomial(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars = save_pars(all = TRUE))

m4n <- brm(bf(amendments ~ vp + words + groups + offset(log(age)), 
              shape ~ vp + groups), 
           data=cndat, 
           family=negbinomial(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars = save_pars(all = TRUE))

m5n <- brm(bf(amendments ~ vp + words + civicAct + offset(log(age)), 
              shape ~ vp + civicAct), 
           data=cndat, 
           family=negbinomial(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars = save_pars(all = TRUE))

m6n <- brm(bf(amendments ~ vp  + words + confid_index + offset(log(age)), 
              shape ~ vp + confid_index), 
           data=cndat, 
           family=negbinomial(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars = save_pars(all = TRUE))

## Table A4
s1 <- summary(m1n)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Gov. Conf.", p=unname(c(sigf(m1n)))) %>% arrange(obs) %>% select(-obs)
s2 <- summary(m2n)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Party Conf.", p=unname(c(sigf(m2n)))) %>% arrange(obs) %>% select(-obs)
s3 <- summary(m3n)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Court Conf.", p=unname(c(sigf(m3n)))) %>% arrange(obs) %>% select(-obs)
s4 <- summary(m4n)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Groups", p=unname(c(sigf(m4n)))) %>% arrange(obs) %>% select(-obs)
s5 <- summary(m5n)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Civic Act", p=unname(c(sigf(m5n)))) %>% arrange(obs) %>% select(-obs)
s6 <- summary(m6n)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Pol. Trust Idx", p=unname(c(sigf(m6n)))) %>% arrange(obs) %>% select(-obs)

all_s <- bind_rows(s1, s2, s3, s4, s5, s6) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "vbl", "p")) %>% 
  mutate(term = gsub("govConfid|courtConfid|partyConfid|civicAct|groups|confid_index", "Social Capital", term), 
         sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, vbl, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") %>% 
  pivot_wider(names_from="vbl", values_from="vals") %>% 
  select(-param)

all_s$term[seq(2,14,by=2)] <- ""

all_s <-all_s %>% 
  mutate(term = gsub("vp", "Constitutional Rigidity", term), 
         term = gsub("words", "log(# Words)", term), 
         term = gsub("shape_", "Shape: ", term))



l1 <- loo(m1n, moment_match=TRUE)
l2 <- loo(m2n, moment_match=TRUE)
l3 <- loo(m3n, moment_match=TRUE)
l4 <- loo(m4n, moment_match=TRUE)
l5 <- loo(m5n, moment_match=TRUE)
l6 <- loo(m6n, moment_match=TRUE)

l_ests <- rbind(l1$estimates[3,], 
                l2$estimates[3,],
                l3$estimates[3,],
                l4$estimates[3,],
                l5$estimates[3,],
                l6$estimates[3,])


loos <- c("LOO IC (SE)", apply(l_ests, 1, function(x)sprintf("%.1f (%.1f)", x[1], x[2])))
names(loos) <- names(all_s)
all_s <- rbind(all_s, loos)
r <- c("R(y, y-hat)", sprintf("%.2f", c(cor(fitted(m1n)[,1], cndat$amendments),
                                        cor(fitted(m2n)[,1], cndat$amendments),
                                        cor(fitted(m3n)[,1], cndat$amendments),
                                        cor(fitted(m4n)[,1], cndat$amendments),
                                        cor(fitted(m5n)[,1], cndat$amendments),
                                        cor(fitted(m6n)[,1], cndat$amendments))))
all_s <- rbind(all_s, r)
all_s

## Calculate First Differences
f1a <- fd(m1n, "vp", cndat, level=0.8) %>% mutate(model = "Confidence in\n Government", 
                                                  variable = "Constitutional Rigidity")
f1b <- fd(m1n, "govConfid", cndat, level=0.8) %>% mutate(model = "Confidence in\n Government", 
                                                         variable = "Social Capital")

f2a <- fd(m2n, "vp", cndat, level=0.8)%>% mutate(model = "Confidence in\n Parties", 
                                                 variable = "Constitutional Rigidity")
f2b <- fd(m2n, "partyConfid", cndat, level=0.8)%>% mutate(model = "Confidence in\n Parties", 
                                                          variable = "Social Capital")

f3a <- fd(m3n, "vp", cndat, level=0.8)%>% mutate(model = "Confidence in\n Courts", 
                                                 variable = "Constitutional Rigidity")
f3b <- fd(m3n, "courtConfid", cndat, level=0.8)%>% mutate(model = "Confidence in\n Courts", 
                                                          variable = "Social Capital")

f4a <- fd(m4n, "vp", cndat, level=0.8)%>% mutate(model = "Group Membership", 
                                                 variable = "Constitutional Rigidity")
f4b <- fd(m4n, "groups", cndat, level=0.8)%>% mutate(model = "Group Membership", 
                                                     variable = "Social Capital")

f5a <- fd(m5n, "vp", cndat, level=0.8)%>% mutate(model = "Civic Activism", 
                                                 variable = "Constitutional Rigidity")
f5b <- fd(m5n, "civicAct", cndat, level=0.8)%>% mutate(model = "Civic Activism", 
                                                       variable = "Social Capital")

f6a <- fd(m6n, "vp", cndat, level=0.8) %>% mutate(model = "Pol. Trust\n Index", 
                                                  variable = "Constitutional Rigidity")
f6b <- fd(m6n, "confid_index", cndat, level=0.8)%>% mutate(model = "Pol. Trust\n Index", 
                                                           variable = "Social Capital")

fd_all <- bind_rows(f1a, f1b, f2a, f2b, f3a, f3b, 
                    f4a, f4b, f5a, f5b, f6a, f6b) %>%
  mutate(credible = factor(ifelse(pval > .9, 2, 1), levels=1:2, labels=c("No", "Yes")))

## Figure A1
fd_all %>% 
  mutate(vbl = factor(vbl, 
                      levels=c("mean", "sd"), 
                      labels=c("Mean", "Std. Deviation"))) %>%
  group_by(model) %>% 
  mutate(so = difference[cur_data()$variable == "Social Capital" & 
                           cur_data()$vbl == "Mean"], 
         model = case_when(model == "Group Membership" ~ "Group\nMembership", 
                           TRUE ~ model)) %>% 
  ungroup %>% 
  mutate(model = reorder(as.factor(model), so, mean)) %>% 
  ggplot(aes(x=difference, y=model, xmin = lwr, xmax=upr, colour=credible)) + 
  geom_pointrange() + 
  geom_vline(xintercept = 0, linetype=3) + 
  facet_grid(variable~vbl) + 
  theme_bw() + 
  theme(panel.grid = element_blank(), 
        legend.position = "top", 
        plot.caption.position = "plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  labs(x="First Difference in Predicted Number of Amendments\nMean and Standard Deviation", 
       y="", 
       colour= "Credible\nDifference (90%)", 
       caption="See Table A4 (models 1-6) for model results ")

ggsave("figures/fd1.png", height=6.5, width=8, units="in", dpi=1000)

## Figure 1
fd_all %>% 
  slice_tail(n=12)%>% 
  mutate(variable2 = case_when(variable == "Social Capital" ~ model, 
                               TRUE ~ variable), 
         variable2 = case_when(variable2 == "Pol. Trust\n Index" ~ "Pol. Trust Index", 
                               TRUE ~ variable2), 
         variable2 = factor(variable2, levels=c("Constitutional Rigidity", "Group Membership", "Civic Activism", "Pol. Trust Index")),
         vbl = factor(vbl, levels=c("mean", "sd"), 
                      labels=c("Mean", "Std. Deviation")), 
         model2 = factor(model, levels=c("Group Membership", "Civic Activism", "Pol. Trust\n Index"), 
                         labels=c("Model 1", "Model 2", "Model 3"))) %>% 
  ggplot(aes(x=difference, y=variable2, xmin = lwr, xmax=upr)) + 
  geom_pointrange() + 
  geom_vline(xintercept = 0, linetype=3) + 
  facet_grid(model2~vbl, scales="free_y") + 
  theme_bw() + 
  theme(panel.grid = element_blank(), 
        legend.position = "top", 
        plot.caption.position = "plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  labs(x="First Difference in Predicted Amendments:\nMean and Standard Deviation", 
       y="", 
       caption="Note: See Table A4 (models 4-6) for model results") 
ggsave("figures/fig1.png", height=8, width=8, units="in", dpi=1000)

## Heteroskedastic Linear Models

m1l <- brm(bf(I(amendments/age) ~ vp + words + govConfid, 
              sigma ~ vp + govConfid), 
           data=cndat, 
           family=gaussian(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207,
           save_pars = save_pars(all=TRUE))


m2l <- brm(bf(I(amendments/age) ~ vp + words + partyConfid, 
              sigma ~ vp + partyConfid), 
           data=cndat, 
           family=gaussian(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207,
           save_pars = save_pars(all=TRUE))
m3l <- brm(bf(I(amendments/age) ~ vp + words + courtConfid, 
              sigma ~ vp + courtConfid), 
           data=cndat, 
           family=gaussian(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207,
           save_pars = save_pars(all=TRUE))

m4l <- brm(bf(I(amendments/age) ~ vp + words + groups, 
              sigma ~ vp + groups), 
           data=cndat, 
           family=gaussian(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207,
           save_pars = save_pars(all=TRUE))

m5l <- brm(bf(I(amendments/age) ~ vp + words + civicAct, 
              sigma ~ vp + civicAct), 
           data=cndat, 
           family=gaussian(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207,
           save_pars = save_pars(all=TRUE))

m6l <- brm(bf(I(amendments/age) ~ vp  + words + confid_index, 
              sigma ~ vp + confid_index), 
           data=cndat, 
           family=gaussian(),
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207,
           save_pars = save_pars(all=TRUE))

## Table A7
s1 <- summary(m1l)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Gov. Conf.", p=unname(c(sigf(m1l)))) %>% arrange(obs) %>% select(-obs)
s2 <- summary(m2l)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Party Conf.", p=unname(c(sigf(m2l)))) %>% arrange(obs) %>% select(-obs)
s3 <- summary(m3l)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Court Conf.", p=unname(c(sigf(m3l)))) %>% arrange(obs) %>% select(-obs)
s4 <- summary(m4l)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Groups", p=unname(c(sigf(m4l)))) %>% arrange(obs) %>% select(-obs)
s5 <- summary(m5l)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Civic Acts", p=unname(c(sigf(m5l)))) %>% arrange(obs) %>% select(-obs)
s6 <- summary(m6l)$fixed %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% mutate(obs = c(1,5,2,3,4,6,7), vbl = "Conf. Index", p=unname(c(sigf(m6l)))) %>% arrange(obs) %>% select(-obs)



all_l <- bind_rows(s1, s2, s3, s4, s5, s6) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "vbl", "p")) %>% 
  mutate(term = gsub("govConfid|courtConfid|partyConfid|civicAct|groups|confid_index", "Social Capital", term), 
         sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, vbl, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") %>% 
  pivot_wider(names_from="vbl", values_from="vals") %>% 
  select(-param)

all_l$term[seq(2,14,by=2)] <- ""

all_l <-all_l %>% 
  mutate(term = gsub("vp", "Constitutional Rigidity", term), 
         term = gsub("words", "log(# Words)", term), 
         term = gsub("sigma_", "Sigma: ", term))


l1 <- loo(m1l, moment_match=TRUE)
l2 <- loo(m2l, moment_match=TRUE)
l3 <- loo(m3l, moment_match=TRUE)
l4 <- loo(m4l, moment_match=TRUE)
l5 <- loo(m5l, moment_match=TRUE)
l6 <- loo(m6l, moment_match=TRUE)

l_ests <- rbind(l1$estimates[3,], 
                l2$estimates[3,],
                l3$estimates[3,],
                l4$estimates[3,],
                l5$estimates[3,],
                l6$estimates[3,])


loos <- c("LOO IC (SE)", apply(l_ests, 1, function(x)sprintf("%.1f (%.1f)", x[1], x[2])))
names(loos) <- names(all_l)
all_l <- rbind(all_l, loos)
r <- c("R(y, y-hat)", sprintf("%.2f", c(cor(fitted(m1l)[,1], cndat$amendments/cndat$age),
                                        cor(fitted(m2l)[,1], cndat$amendments/cndat$age),
                                        cor(fitted(m3l)[,1], cndat$amendments/cndat$age),
                                        cor(fitted(m4l)[,1], cndat$amendments/cndat$age),
                                        cor(fitted(m5l)[,1], cndat$amendments/cndat$age),
                                        cor(fitted(m6l)[,1], cndat$amendments/cndat$age))))
all_l <- rbind(all_l, r)
all_l

## Figure A2
## Calculate first differences

f1a <- fd(m1l, "vp", cndat, dist="normal", level=0.8) %>% mutate(model = "Confidence in\n Government", 
                                                                 variable = "Constitutional Rigidity")
f1b <- fd(m1l, "govConfid", cndat, dist="normal", level=0.8) %>% mutate(model = "Confidence in\n Government", 
                                                                        variable = "Social Capital")

f2a <- fd(m2l, "vp", cndat, dist="normal", level=0.8)%>% mutate(model = "Confidence in\n Parties", 
                                                                variable = "Constitutional Rigidity")
f2b <- fd(m2l, "partyConfid", cndat, dist="normal", level=0.8)%>% mutate(model = "Confidence in\n Parties", 
                                                                         variable = "Social Capital")

f3a <- fd(m3l, "vp", cndat, dist="normal", level=0.8)%>% mutate(model = "Confidence in\n Courts", 
                                                                variable = "Constitutional Rigidity")
f3b <- fd(m3l, "courtConfid", cndat, dist="normal", level=0.8)%>% mutate(model = "Confidence in\n Courts", 
                                                                         variable = "Social Capital")

f4a <- fd(m4l, "vp", cndat, dist="normal", level=0.8)%>% mutate(model = "Group Membership", 
                                                                variable = "Constitutional Rigidity")
f4b <- fd(m4l, "groups", cndat, dist="normal", level=0.8)%>% mutate(model = "Group Membership", 
                                                                    variable = "Social Capital")

f5a <- fd(m5l, "vp", cndat, dist="normal", level=0.8)%>% mutate(model = "Civic Activism", 
                                                                variable = "Constitutional Rigidity")
f5b <- fd(m5l, "civicAct", cndat, dist="normal", level=0.8)%>% mutate(model = "Civic Activism", 
                                                                      variable = "Social Capital")

f6a <- fd(m6l, "vp", cndat, dist="normal", level=0.8) %>% mutate(model = "Pol. Trust\n Index", 
                                                                 variable = "Constitutional Rigidity")
f6b <- fd(m6l, "confid_index", cndat, dist="normal", level=0.8)%>% mutate(model = "Pol. Trust\n Index", 
                                                                          variable = "Social Capital")

fd_all <- bind_rows(f1a, f1b, f2a, f2b, f3a, f3b, 
                    f4a, f4b, f5a, f5b, f6a, f6b) %>%
  mutate(credible = factor(ifelse(pval > .9, 2, 1), levels=1:2, labels=c("No", "Yes")))

## Make Figure
fd_all %>% 
  mutate(vbl = factor(vbl, 
                      levels=c("mean", "sd"), 
                      labels=c("Expectation", "Std. Deviation"))) %>%
  group_by(model) %>% 
  mutate(so = difference[cur_data()$variable == "Social Capital" & 
                           cur_data()$vbl == "Expectation"], 
         model = case_when(model == "Group Membership" ~ "Group\nMembership", 
                           TRUE ~ model)) %>% 
  ungroup %>% 
  mutate(model = reorder(as.factor(model), so, mean)) %>% 
  ggplot(aes(x=difference, y=model, xmin = lwr, xmax=upr, colour=credible)) + 
  geom_pointrange() + 
  geom_vline(xintercept = 0, linetype=3) + 
  facet_grid(variable~vbl) + 
  theme_bw() + 
  theme(panel.grid = element_blank(), 
        legend.position = "top", 
        plot.caption.position = "plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  labs(x="First Difference in Predicted Amendment Rate\nMean and Standard Deviation", 
       y="", 
       colour= "Credible\nDifference (90%)", 
       caption="See Table A5 (models 1-6) for model results")
ggsave("figures/fd2.png", height=6.5, width=8, units="in", dpi=1000)

